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Abstract 

In this paper, we propose an efficient Monte Carlo implementation of non-linear 
FBSDEs as a system of interacting particles inspired by the ideas of branching diffu- 
sion method. It will be particularly useful to investigate large and complex systems, 
and hence it is a good complement of our previous work presenting an analytical per- 
turbation procedure for generic non-linear FBSDEs. There appear multiple species of 
particles, where the first one follows the diffusion of the original underlying state, and 
the others the Malliavin derivatives with a grading structure. The number of branch- 
ing points are capped by the order of perturbation, which is expected to make the 
scheme less numerically intensive. The proposed method can be applied to semi-linear 
problems, such as American and Bcrmudan options, Credit Value Adjustment (CVA), 
and even fully non-linear issues, such as the optimal portfolio problems in incomplete 
and/or constrained markets, feedbacks from large investors, and also the analysis of 
various risk measures. 
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1 Introduction 



The forward backward stochastic differential equations (FBSDEs) were first introduced by 
Bismut (1973) [T], and then later extended by Pardoux and Peng (1990) [27] for general 
non-linear cases. They were found particularly relevant for optimal portfolio and indiffer- 
ence pricing issues in incomplete and/or constrained markets. Their financial applications 
are discussed in details in, for example, El Karoui, Peng and Quenez (1997a) [10J, Ma and 
Yong (2000) [23] and a recent book edited by Carmona (2009) [4]. The importance of 
FBSDEs will increase in coming years even among practitioners where the new financial 
regulations will put significant constraints on available assets and trading strategies. 

In the recent paper, Fujii & Takahashi (2011) [12] proposed a new perturbative so- 
lution technique for generic non-linear FBSDEs. It was shown that a non-linear FBSDE 
can be decomposed into a series of linear and decoupled FBSDEs by treating a non-linear 
driver and feedback terms as perturbations to the corresponding decoupled free system. 
In particular, it allows analytical explicit expressions for the backward components with 
the help of the asymptotic expansion technique (See, for example [281 ED HEB [29].). A 
backward component of the diffusion part was shown to be obtained by directly consider- 
ing dynamics of the stochastic flow, which denotes a Malliavin derivative of the underling 
state process, or simply applying Ito formula to the result of the other part. In Fujii & 
Takahashi (2012) [13], the method was applied to a quadratic-growth FBSDE appearing 
in an incomplete financial market with stochastic volatility. Explicit expressions for both 
of the backward components were obtained up to the third order of the volatility of volatil- 
ity. The comparison to the exact solution with Cole-Hopf transformation demonstrated 
effectiveness of the perturbative expansion. 

Notice the fact that one can already apply standard Monte Carlo simulation to the 
results obtained in each order of the perturbative expansion in [12J . However, due to its 
convoluted nature, it contains multi-dimensional time integrations of expectation values 
which make the naive applications too time consuming, particularly for the evaluation 
of higher order perturbation terms. To handle this problem, we applied the idea of par- 
ticle representation used in branching diffusion models, such as in McKean (1975) [26J. 
There, the convoluted expectation is compressed into a single standard expectation by 
introducing an intensity of the particle interaction. McKean [26] applied the method to 
solve a particular type of semi-linear PDE, where a single particle splits into two at each 
interaction time and creates a cascade of the identical particles. Note that, our method 
is not directly related to McKean [26] since the interested system is already decomposed 
into a set of linear problems, although we have used the similar particle representation to 
avoid nested simulations. 

The analysis of branching Markov process and related problems in semi-linear PDEs 
has a long history. Some of the well-known works are Fujita (1966) [15] . Ikeda, Naga- 
sawa & Watanabe (1965,1966,1968) [UJ Mi EH] , Ikeda et.al. (1996,1997) [20] and Naga- 
sawa & Sirao (1969) [25]. As for a recent work, in particular, Chakraborty & Lopez- 
Mimbela (2008) used particle representation where the number of offspring at each inter- 
action point is randomly drawn by some probability distribution, which can be finitely 
many or infinite. The authors used the branching particle representation to study the 
existence of global solutions for semi-linear PDEs with a non-linear driver given by a 

x The same branching representation is already seen in [19] . for example. 
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generic polynomial function □• Recently, Henry-Labordere (2012) [16] introduced a par- 
ticle representation to study the semi-linear problems in finance. He called it marked 
branching diffusion and has discussed its application to efficiently calculate CVA (credit 
value adjustment) in one-shot Monte Carlo simulation. He also referred to its application 
to other semi-linear problems, such as American options, as well as its possible extension 
to truly non- linear problems by using Malliavin derivatives. 

In the current paper, we combine the idea of particle representation and the per- 
turbation technique developed in the previous work [12J. We provide a straightforward 
simulation scheme to solve fully-nonlinear decoupled as well as coupled FBSDEs at each 
order of perturbative approximation. In contrast to the direct application of branching 
diffusion method, the number of branching points are capped by the order of perturbative 
expansion, which is due to the linearity of the decomposed FBSDE system. This property 
is expected to make Monte Carlo simulation less numerically intensive. Our method can be 
applied to semi-linear problems, such as American and Bermudan options El, Credit Value 
Adjustment (CVA) as special examples. It can be also applied to fully non-linear (and fully 
coupled) issues, such as the optimal portfolio problems in incomplete and/or constrained 
market, analysis for various risk measures as well as for the feedbacks from so-called large 
investors. Concrete applications of the new method will be published separately |14j . 

2 Setup 

We first consider generic decoupled non-linear FBSDEs. Let us use the same setup assumed 
in the work [12]. The probability space is taken as (0, J 7 , P) and T € (0, oo) denotes some 
fixed time horizon. Wt = (Wt-> " ' ■> ^T)*> < i < T is Revalued Brownian motion defined 
on (O, J 7 , P), and (^i){o<t<T} stands for P-augmented natural filtration generated by the 
Brownian motion. 

We consider the following forward-backward stochastic differential equation (FBSDE) 

dV s = -f{X s ,V s ,Z s )ds + Z s -dW s (2.1) 
V T = ¥(Xr) (2.2) 

where V takes the value in R, and Xt € M. d is assumed to follow a generic Markovian 
forward SDE 

dX s = j (X s )ds + j(X s ) ■ dW s . (2.3) 

Here, we absorbed an explicit dependence on time to X by allowing some of its components 
can be a time itself. \P(Xy) denotes the terminal payoff where ^(x) is a deterministic 
function of x. The following approximation procedures can be applied in the same way 
also in the presence of coupon payments. Z and 7 take values in R r and M. dxr respectively, 
and " ■" in front of the dW represents the summation for the components of r-dimensional 

2 For recent developments and reviews of the particle methods, see for examples 8, 9 . There exist a 
significant amount of works related to branching diffusion in 1960's and 70's. There are also a vast range 
of new applications and enhancements in biology, such as gene mutation and population growth problems, 
as well as in engineering issues. We have not yet obtained the whole picture of research history related to 
branching diffusion and are welcoming information from those familiar with the topic. 

3 A BSDE formulation for an American option was shown in El Karoui etal. (1997b) [TT], which was 
recently studied by Labart & Lelong (2011) [22] based on regression based Monte Carlo simulation. 
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Brownian motion. Throughout this paper, we are going to assume that the appropriate 
regularity conditions are satisfied for the necessary treatments. 

Let us fix the initial time as t. We denote the Malliavin derivative of X u (u > t) at 
time t as 

V t X u € R rxd (2.4) 

Its dynamics in terms of the future time u is specified by the well-known stochastic flow: 

d(Y t>u )i = d kl i{X u ){Y tu ))du + d kl i{X u ){Y tu ))dW* 

(Yt,t)} = % (2.5) 

where d k denotes the differential with respect to the k-th component of X, and 5j denotes 
Kronecker delta. Here, i and j run through {1, • • • ,d} and {1, ■ ■ ■ , r} for a. Throughout 
the paper, we adopt Einstein notation which assumes the summation of all the paired 
indexes. Using the known chain rule of Malliavin derivative, one sees 

(v t xi)= d k y Q (x s )(v t x*)ds + / d k Y(x s )(v t x*)-dw s + f(x t ) (2.6) 

and hence it satisfies 

(V t Xi) a = (Yt^jUXt) = (Y tiUl (X t ))i (2.7) 
where "a" is the index of r-dimensional Brownian motion. 



3 Expansion into a series of Linear FBSDE System 

Following the perturbative method proposed in [12], let us introduce the perturbation 
parameter e and then write the equation as 

J dV s {e) = -ef(X s ,V s {e \z { s e) )ds + Z ( s e) -dW s 

{ v¥ = *{x T ) {3A) 

where e = 1 corresponds to the original model 0. We suppose that the solution can be 
expanded in a power series of e: 

V ^ = f + e y i (1) + e V/ 2) + e V< 3) + ... (3.2) 

Zf = zf ) + eZ[ l) +e 2 z[ 2) +e*Z ( ? ) + --- (3.3) 

If the non-linearity is sub-dominant, one can expect to obtain reasonable approximation 
of the original system by putting e = 1 at the end of calculation. 

The dynamics of each pair (V®, Z^) can be easily derived as follows: 
Zero-th order 



4 It is possible to extract the linear term from the driver and treat separately. Here, we simply leave it 
in a driver, or work in a "discounted" base to remove linear term in V. 
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First order 



Second order 



dV s {1) = -f(X s , Vp,Zp)ds + zi l) ■ dW s 
V® = 



dvP = - {v^l + (Z^)^r} fix., vP,z^)ds + Z® ■ dW s 



i4 2) = o 



Third order 



)\2 d 2 xU (l)ya(l)j2 



s dvdz a 



= 



+\z a s {1) z b s {1) ^}f{Xs, vF>,zF)ds + z® . dWs 



One can continue to an arbitrary higher order in the same way. 



(3.5) 



(3.6) 



(3.7) 



Note that the higher order backward components (V^ n \ Z( n )){ n>1 } are always outside 
of the non-linear functions. This property arises naturally due to the very nature of 
perturbation. As we shall see, this is crucial to suppress the number of particles in the 
numerical simulation. 

4 Interacting Particle Interpretation 

Let us fix the initial time t and set Xt = xt- 



4.1 e-Oth Order 

For the zeroth order, it is easy to see 



(0) 



Z 



o(0) 



*(X T ) 



d^(X T )(VlX l T ] 



Ft 



d^{X T ){Y tTl {X t )t 



F 



(4.1) 



(4.2) 



It is clear that they can be evaluated by standard Monte Carlo simulation. However, 
for their use in higher order approximation, it is crucial to obtain explicit approximate 
expressions for these two quantities. As proposed in |12j . we use asymptotic expansion 
technique [28j EU [30j [29] for this purpose. When ^ is a smooth function, it is quite 
straightforward. Even if ^ is not a smooth function, such as an option payoff, one can ob- 
tain explicit expressions of (V(°), Z<°)) in terms of X t , too. This is because, one can derive 
an approximate joint transition density of general diffusion processes by the asymptotic 
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expansion In the following, let us suppose that we have obtained the solutions up to 
given order of asymptotic expansion, and write each of them as a function of xt- 



(o) 



(4.3) 



4.2 e-lst Order 

Since the BSDE is linear, we can integrate as before. Here, let us first consider the 



evaluation of V+ 



(i) 



(i) 



E 



E 



f(x u ,v(°\zW) 



f[X u ,v^(X u ),z^(X u ) 



T t 



(0), 



du 



T t 



du 



(4.4) 



Although it is possible to carry out standard Monte Carlo simulation for every time u € 
(t, T) and integrate to obtain the V t , the time integration becomes numerically quite 
heavy. In fact, it will soon become infeasible for e higher order terms that include multi- 
dimensional integration of time. We now introduce particle interpretation by McKean 
developed for the study of semilinear PDEs: 

Proposition 1 The in b4-4\) can be equivalently expressed as 



(i) 



l{r>t}IE 



i {T<T} ft(x T M°\x T ),z^(x T ) 



T 



(4.5) 



Here r is the time of interaction which is drawn independently from Poisson distribution 
with an arbitrary deterministic positive intensity process Xt- It can be a positive constant 
for the simplest case. ^ f is defined as 



f t (x,v ( -°\x),z^(x)) 



(0), 



f(x,v^{x),z^(x)) 



(0), 



(4.6) 



Proof: Define the new process for (s > t): 



then its dynamics is given by 



Since we have V$ = V^\ one can easily see the following relation holds: 



(4.7) 



e f t *Kdu {\ sV W ds _ f(X s ,v (0 \X s ),z^\X s ))ds + ZP ■ dW s } 
XsV^ds - \J t (X s ,v(°\X s ),z(°\X s ))ds + eft x " du zV . dWs . (4i 



(i) 



E 



-K x ° ds \ u f t (X u ,vW(X u ),z^(X u )) 



T t 



du 



5 We intend to use the result of asymptotic expansion only for higher order approximations. 
6 It is not difficult to make it a stochastic process. 



(4.9) 
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It is clear for those familiar with credit risk modeling [21 [3], it is nothing but the present 
value of default payment where the default intensity is A with the default payoff at s (> t) 
as f t (X s ,v(°\X s ),z(°\X s )). Thus, it is clear that flO) is equivalent to gSJ). ■ 



Now, let us consider the martingale component Z^ l \ It can be expressed as 



.(1) 



E 



V t f(x u ,v(°\X u ),z(°\X u ] 



du 



(4.10) 



We perform the similar transformation for Z^ 1 ' to make it easier to interpret in the inter- 
acting particle model. Firstly, let us observe that the dynamics of Malliavin derivative of 
follows 

d{V t VP) = -(V t Xl){d l + dM°\X s )d v + d l z a ^\X s )d z a}f(X s ,v^\X s ),z^\X s ))ds 

+ (V t ZP) ■ dW s (4.11) 

V t V t W = zf ] (4.12) 
For lighten the notation, let us introduce a derivative operator 

Vi(x,v(°\zW) = di + div^{x)d v + d iZ a{Q \x)d z a (4.13) 

and also 

f(x,v^,z^) = f(x,v^(x),z^(x)) (4.14) 
Now, we can write Eq. (I4.12j) as 

d(V t vP) = -(V t Xi)Vi(X s ,v(°\zW)f(X s ,v(°\zW)ds + (V t zP) ■ dW s 
Define, for (s > t), 

V t V s W = eft x ^ du (V t V s m ) (4.15) 
then its dynamics can be written as 

d(V t V s {1) ) = eft x - du {\ s (V t VP)ds - (VtXi)Vi(X 8 M \z®MX.M°\* {0) )da 
+V t Z® ■ dW s ] 

= X s {V t vP)ds - UVtXDViiX^v^^^ftiX^v^^^ds 

+e s t s ^ du (v t z^)-dw s 

We have 



v t v t {1) = z[ 1] 



and hence 



E 



e-^^ ds \ s {V t Xi)V i {X u ^),z^)f t {X u M°\z^) 



Ft 



(4.16) 
(4.17) 
(4.18) 



Thus, following the same argument of the proposition [H we can conclude: 
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Proposition 2 in is equivalently expressed as 



70.(1) 



l {r>t} 



E 



l { r < T } (Yt,r7(X t )tV l (X T y°\z^)f t (X T ,V^,Z^) 



(4.19) 



where the definitions of random time r and the intensity process A are the same as those 
in proposition [0 

As we shall see later, interpreting (X, Y) as a pair of particles allows an efficient Monte 
Carlo implementation. For the evaluation of for example, one can consider it as an 
system of two particles (X, Y), which have the intensity A of the interaction that produces 



(Y tT rtX t ))iV i (X T ,v(°\zW)fc(X T ,v( \zM) 



(4.20) 



at the interaction point and annihilate altogether. For the interpretation is much 

simpler. A single particle X with the decay rate of A leaves / at its decay point and 
vanishes. 



4.3 e-2nd Order 

For the e-2nd order, one can observe that 

i-T 



(2) 



(2) 



E 



E 



v^d v + z^d za )f(x u M°\z^) 



du 



v t {(y^d v + z^d za )f(x u M°\^)} 



T 



du 



(4.21) 
(4.22) 



solve the BSDE (|3.6p . Its particle interpretation is available by the similar transformation. 
Firstly, for (s > t), let us define 



Vff =e i:^uy(2) 

with some appropriate intensity process A. Then it follows 



(4.23) 



dV t 



(2) 



X s V t fds - X s (V s ^d v + Z^d za )f t (X s , v<-°\zW)d6 

+e ft *udu z (2) . dw 



" (2) (2) 

Observing that V t t =V t , one can confirm that 



V, 



(2) 



l{n>*} E 



l {T1<T} (v^d v + Z^d za )f t (X T1 M°\z^) 



T 



(4.24) 



(4.25) 



where t\ is the random interaction time with intensity A. Now, using the tower property 
of conditional expectations, one can conclude that 



(2) 

Proposition 3 V t in h4-21\) is equivalently expressed as 



V} 2) = l {Tl>t} E 



<T2 <T} \pv /t,Tl ) fn ,T2 
+ 1 {ri>t} E 1 {ri<r 2 <T} , 92 a /t,T 1 (Y Ti ,t 2 1tA ^i,T 2 f 



T 



(4.26) 
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where we have defined 

ft,s = ft(XsM Hx a ),z(°\x a )) 

V hS ^V i (X s ,v(°\X s ),z(°\X s )) 

It = l(X t ) (4.27) 

and T\ and T2 are the two interaction times randomly drawn with intensity A. 

A particle interpretation for the first term is quite simple. A particle X starts at t 
follows the diffusion (j2.3j) with (self) interaction intensity A. For the first interaction time 
T\, it yields d v ft )T1 and at the 2nd interaction time T2 it yields f T1 ,T 2 an d decays away. The 
expectation value can be evaluated by preparing a large number of particles X starting 
from the same point and obeying the same physical law but spend independent lives. For 
the second term, the interpretation is more interesting. A particle X starts at time t and 
follows the diffusion (|2.3h with interaction intensity A. At the first interaction time t\, it 
yields d z af t>Tl and at the same time bears a new particle Y. After n, the two particles 
(X,Y) follow the diffusions (12. 3p and (12. 5p . respectively. They have interaction intensity 
A, and at the second interaction point T2 they yield (Y TljT2 j(X T2 ))S/ T2 f T1:T2 and annihilate 
altogether. As in the first example, the expectation can be calculated by preparing a large 
number of particle X at the same starting point. 



Remark: Note that, if we simply use Eqs. {4-4\ \4-®l ) an d the tower property, we have 



to handle a two-dimensional time integration. It makes naive implementation of Monte 
Carlo simulation numerically too heavy. In our particle interpretation, this problem is 
solved by introducing random interaction times with some intensity A. One can choose 
appropriate size of intensity that produces enough amount of events for Monte Carlo sim- 
ulation. 

We now consider an interacting particle interpretation of Z^ 2 \ For the evaluation of 
Z( 2 \ we need to define the second order stochastic flow for (t < s < u): 

ox t ox* ox t 

Since we have 

{Y s>u ){ = 51+ [ U (Y s J k (d n l (X v )dv + d a \X v ) ■ dW v ) (4.29) 

J s 

it is easy to see that 

(r t , s ,4, fc = J {v t , s , v ) l jjk {dai{x v )dv + d n l (x v ) ■ dw v ) 

+ [ {Y tiV )f(Y S!V )[{d lm %{X v )dv + d lm Y{X v )-dW v ) (4.30) 



Note that we have Tt tS ^ s = 0, regardless of time s (> t). Using the second order stochastic 
flow, the Malliavin derivative of Y can be written as 

V?(Y s , v )l = (T tjS , v )i k (ji(X t )) a = (r t , s , u7 pQ))L (4-31) 
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(2) 

Consider the process of Malliavin derivative VtVg ■ One can write its dynamics for 
(t < s) as 



d(V t vP) = -[{V t vP)d v + (V t zfV)d z «)f(X s ,v(°\zW)ds 

+V t ZP ■ dW s (4.32) 



VtV. 



(2) 



z 



(2) 



As before, we define 

v t vP = e ft Kdu (v t vP) 

then its dynamics satisfies the following SDE: 



(4.33) 



(4.34) 



d(V t vP) = \ s {V t vP)ds - X s (VtX^lvPViJdJtJ + (Z^)V l , s (d za f t , s ) 



r(2)> 



+ ({V t VP)d v + (V t Z^)d z a)f t Jds + e ft x " du V t zW . dW ; 



(v t v} 2) ) = z® 

Then, the same arguments leads to 



(4.35) 
(4.36) 



(2) 



l{n>i} E 



l^TyiVtX^fv^Wi^dJt^) + (Z<V)Vi, Tl (d z af t>Tl ] 



+l {T1<T} [(V t V T W)d v + {V t Z«U)d za )f\ T1 
using the random interaction time t\. 

Proposition 4 Z\ ' in [4-%ty ^ s equivalently expressed as 



F 



(4.37) 



Z 



o(2) 



l{n>t} E 



1 {ri<T2<T}(^,T 1 7t)a V j,Ti(^/t,T 1 )/ri,T2 



+ l{r 1 >t} E 1 {T 1 <r 2 <T}(^,r 1 7i)aVi,r 1 (<9^/t, T1 )(y T1 , T2 7r 1 )fcV i)T2 / Tlir2 



+ 1 {r 1 >t} E 1 {rt<7s<T}(^«/t,n)( i t,757t)o^i.75/n,75 



■Ft 



+ 1 



+ 1 



+ 1 



{n>t} 



E 



1 {T 1 <T2<T}(dz b ft,T 1 )(7T 1 )l(rt,T 1 ,T27t)) t a^i,T 2 fT 



Ft 



{n>t} 



{n>t} 



E 



L {ti<t 2 <T} ,ri7t)a(^j'7ri )b (^ / ri,T 2 )fc^' 'i,T2 /ri,7 



F 



E 



1 {ti<t 2 <T} ,T 2 7t)a 

(^r 1 ,r 2 7r 1 )feV i)T2 (V 4iT2 / TliT2 ) LF t . (4.38 



where T\ and T2 are sequential interaction times with intensity A. 

Proof: It can be shown straightforwardly by using the tower property of conditional expec- 
tations and commutativity between the indicator functions and the Malliavin derivative 
due to the independence of A. ■ 
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T 2 < T 




Despite the apparent complexity, required numerical procedures for the evaluation of 
is, in fact, quite simple. We provide a Feynman diagram for the particle interpreta- 
tion in Figure [TJ At the first stage, there are two particles of (X.,Yt t .) with initial values 
(xt, {£*•}), which survive until the second interaction time T2 (< T). At the first inter- 
action at n, two additional particles (Y Tlt .,Tt T i,-) are crea ted. Each interaction occurs 
randomly with intensity A. Note that we already know the initial values of the new par- 
ticles regardless of the interaction time, which makes numerical simulations possible to 
carry out. What one has to do is to store the information of t\ and T2 and the values of 
the particles at these times. Then, all the ingredients in expectations can be calculated. 
Simply repeating independent experiments and taking average will give the desired values. 

4.4 e-3rd Order: 

In the similar fashion, we can proceed to higher order. As before, by considering the 
dynamics of 



(4.39) 



one can observe that 



V® = l {T1>t} E 



i {T1<T} (v^d v + z^d za + Uv^fdl 



(4.40) 



It can be written in terms of the fundamental variables simply applying tower property. 
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(3) 

Proposition 5 V t can be expressed as 



(3) 



l {ri>t} 



E 



- L {ti<T2<T 3 } (d v ft 



T2,T3 



+ 1 



+ 1 {n <r 2 <r 3 } /« ,n ) I (^n ,ra 7n )a V i |T2 (9„ f n )T2 ) / 

+ ( , T2 7ri ) a ^ « , T2 ( ^2 6 /n ,72 ) ( ^r~2 , T3 7t2 ) ft V j ( r3 / r2 ) T:j 
ti,t 2 ,t- a 1t\ )j',aV 1,73/72 ,73 
,t 2 7ri )a iPjlvi )b (^"2 ,T3 )jt Vi jT3 / T2 ,T3 
7r 1 )i(^r 2 ,r 3 77,) l 6 V^(V^/ r2)73 )} 

. 2 
1 {ri<T}^(^/t,T 1 ) ( l{ n<7 f ^/n.Tf 



(4.41) 



{n>*} 



E 



P =i 



/ - \P=1 / a \P=2 

+ 1 {n <T} /t, ri ) y 1 { n < r | <t} fn ,r 2 p J 1 { n <t| <T} ,r 2 p 7n )a V i)7 | / n )7 | t 

1*/ .a \p=i 

+ 1 {r 1 <T} ^ (d z * z bft, T1 ) ( !{ri <r 2 p <T} (Xn ,rf 7n )a V i)7 | / n j7 | 



1 {ti <rf <T} (Xi ,r| 7n )& V jjT | / n >r . 







^p=2 


•Ft 



(4.42) 



where the contents within each bracket of p € {1,2} musi 6e calculated according to the 
diffusion processes {X. n , 5^ 1 ,.) p ={i j 2} that follow the identical diffusion laws with the same 
initial values, but are independent with each other. {rj}j>i are sequential random times 
of interactions drawn with intensity A. {i~2}p=i,2 should be drawn independently. 

Note that, we have introduced two sets of particles labeled by p G {1, 2} that follow the 
same physical laws but perfectly independent with each other to eliminate ri-conditional 
expectations. In this way, one can avoid the use of nested Monte Carlo simulation. In 
Figures [2] and El we have provided two Feynman diagrams, one for the first half, and 

(3) 

the other for the second half of the expression of . In simulations, one has to store 
the interaction times and all the relevant particles values at those points to evaluate the 
expectations. 



r 3 <T 




(3) 

Figure 2: A particle interpretation for the first half of V t . 



12 




< T 



Figure 3: A particle interpretation for the second half of . 
4.5 and e-higher order terms 

The valuation procedures for Z^ are almost the same as that of Z^ 2 \ but we need to 
introduce a new type of particle corresponding to (-£^-T StUyVy .). As easily guessed from 
the previous examples, we need to add one new particle corresponding to a higher order 
stochastic flow to complete the particle picture at every time when we proceed a e-higher 
order approximation (of martingale component Z). A remarkable fact is that all the 
initial conditions of the new particles created at random times are known beforehand 
thanks to the characteristics of the Malliavin derivatives. This feature makes one can 
perform numerical simulations that describe full history of the evolution of particles. 



5 Extension to Fully- Coupled Cases 

We now consider the situation where the underlying state process X also gets the feedbacks 
from the backward components. By making use of the perturbative technique in PDE 
framework [12] . we shall show that the same strategy in the previous sections works well 
also in this seemingly much more complicated situation. 
The dynamics of whole system is given by 



dV t = -f{t, X t , V t , Z t )dt + Z t ■ dW t 
V T = *(Xr) 

dX t = 7o (t, X t , V t , Z t )dt + 7 (t, X t , V t , Z t ) ■ dW t 
X = x 



(5.1) 



where we have distinguished time arguments from X to make PDE generator a familiar 
form. As before, we assume that V, Z,X take value in M,R r and M. d respectively, and W 
denotes a r-dimensional Brownian motion. Following the idea of four-step scheme [23] . we 
postulate that Vt is given by some appropriate function of t and X, v(t, X). Then it needs 
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to satisfy the relevant PDE: 

d t v(t,x) + ldiv(t,x)^Q(t,x,v(t,x),z(t,x)) + \dijv(t,x){^ 1 ■ -f j )(t,x,v(t,x),z(t,x))^ 

+f(t, x, v(t, x),z(t, x)) =0 ^ ^ 

z(t, x) = div(t, x)7*(t, x, v(t, x), z(t, x)) 
v(T,x) = *(T,x) 

The above non-linear PDE cannot be solved in general. Therefore, let us introduce 
perturbation parameter e as before, 



dV, 



(s) 



ef{t,X , f\v t {e \z i f ) )dt + Zl t >-dW i 



dX 



(«0 



(<--) 



r(t,X^) + e f i(t,XrX t ',Zn)dt 



(0 t/( £ ) 



(0 T/( £ ) *-(«0' 



and its corresponding PDE 

^( e )(M) + {d^( e )(t,z^ 

+e/(t,x,t;( £ ),z( e )) = 
z^(t,x) = a^( e )(t,x)7 i (t,x,t;( e ),z( e )) 

w ( £ )(r,x) = *(jc) 

Here, we have extracted the terms free from feedback effects from X's dynamics 0: 

J 7o(t,x,t>^,z( e )) = r(t,x) + €fi(t,x,v^(t,x),z^(t,x)) 
\ 7(t,x,v (e) ,z (e) ) = cr(t,x) + ei](t,x,v^ (t,x), z^ (t,x)) 



(5.3) 



(5.4) 



We suppose that the solution of the above PDE can be expanded perturbatively in 
terms of e as 



v ( e \t,x) = v^(t,x) + evW(t,x) + eV 2 )(i,x) H 

z^) (t, x) = 2fC°) (f , x) + ez (1) (i, x) + e 2 z (2) (f, x) H 



(5.5) 
(5.6) 



As in the previous sections, putting e = 1 is expected to give the approximation of the 
original system as long as the non-linear effects are perturbative. 



5.1 Expansion of non-linear PDE 

Straightforward calculation allows us to expand the original PDE into a series of linear 
parabolic PDEs. See |12] for details. Firstly, let us define the differential operator L: 

C(t,x) = r\t,x)di + ■ ^)(t,x)dij (5.7) 

7 Although this can be done somewhat arbitrarily, it may be natural to set r(t, x) and a(t, x) as the 
expected dynamics of X when all the feedback effects are switched off. 
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which corresponds to the infinitesimal generator of X^°\ ie., the free forward component 



dxf ] = r(t, X t )dt + a(t, X t ) ■ dW t (5.8) 
X^ 0) = x (5.9) 

Using this generator, we can show that the backward components in each order satisfy: 
Zero-th order 

f (d t + £(t,x)y°\t,x) = 
\ v(°\T,x) = V(x) 

and 

,x)a l (t,x) (5-11) 

Higher expansion order (n > 1) 

% + C(t, xfj (t, x) + G< n > (t, x) = 
v( n \T,x) = 



(5.12) 



where the expression of G(") and can be obtained straightforwardly by extracting 
0(e n ) terms from ((531) . 



5.2 Particle Interpretation 

The crucial point in the previous subsection is, because of the perturbation structure in 
(|5.4p . the relevant differential operator always derived from X^ and remains the same 
for all the expansion orders. In addition, since we put a e-factor in front of the non- 
linear terms, contains the backward components with e-order only up to (n — 1). 
Furthermore, it is clear to see that z^ n ' can only contain the backward components of 
{v( m )}| m<n | and {z( m ) }{m<( n _i)}- Therefore, using Feynman-Kac theorem, we see that 
the PDE in (|5.12p is equivalently expressed by 

{ vP = o (5 ' 13) 

where the dynamics of the forward component X^ is already known. Because of the 
very nature of the perturbative expansion, all the (yN^Wjj 

m>i} appear as a power 

(n) 

series and not contained within the non-linear functions. Thus, V t can be solved by the 
same procedures studied in the previous sections, and also the nice properties of explicitly 
capped number of branches and interaction points still hold. 

Note that z[ n ^ is not equal to Z^ that contains additional terms through the feedbacks 
to X. However, it is not difficult to calculate these terms. For example, one can observe: 
1st order (n = 1) 

G {1) (t,x) = f (0) (t,x) +d iV ^(t,x)^°\t,x) +d ijV (° ) (t,x)(a i - V j ^)(t,x) (5.14) 
z (1) (t,x) = d i v il) (t,x)a i (t,x) + d i v W (t,x)i 1 m (t,x) (5.15) 
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2nd order (n = 2) 

G&(t,x) = (vW(t,x)d v + z a W(t,x)d za )f(°\t,x) 

+ 8 iV ^ (t, (t, X) + diV^ ( t , x) («W ( tj x )d v + z a ^ (t, X)d z a)^ (t, x) 



z^(t,x) 



+d ij v^(t,x){a i ■ rf^)(t,x) + ^d ij v^(t,x)(7 ] i ^ ■ rf^)(t,x) 
+d ijV ^ (t, x)a l (t, x) ■ [v {1) (t, x)d v + z a(1) (t, x)d z *)r?W ( tj x ) 
0^ (2) (t, x) + (t, x)rf (0) (t, x) 

+cW 0) (t, x) (t, x)0„ + z a{1) (t, x)d z a)rfW (t, x) 



(5.16) 
(5.17) 



Higher order cases can be obtained similarly. 

Let us now consider the particle method to evaluate the relevant terms. Let us fix 
the initial time as t as before: For the zero-th order, the problem is exactly the same as 
the decoupled case and we can derive easily y(°\t, x) and z^(t, x) as a function of x by 
asymptotic expansion 0. For simplicity, we write as X s , since the underlying process 
does not change. 



1st order 

As for the first order, observe that G^(t,x) is given as an explicit function of x after the 
completion of the zero-th order calculation. Then, vj 1 ^ follows 



dV s {1) = -G^(t,X s )ds + Z^ ] ■ dW, 



(i) 



and hence, by the same arguments, for (s > t), we have a particle representation as 



(5.18) 



V t (l) = l {T>t} E\l {T<T} G ( t 1) (r,X T 



where G^ is defined as 



gS 1} ( S ,X s ) = ^^ du G^(s,X s ) 



(5.19) 



(5.20) 



with some appropriate positive deterministic (or independent) intensity A. For martingale 
component, it is easy to see 



Z 



(i) _ 



Z\ 1] +d i vW{t,x t )rfW{t,x t ) 



(5.21) 



from (|5.15p . Here, the particle representation of Z^ can be derived in the same way as 
in the decoupled case: 



Xi) 



l{r> f } E 



l {T<T} (Yt,T<Tt) t adi&? :) (T,X T 



(5.22) 



As before, this is only to use higher order expansion. For the valuation of the zero-th order itself, one 
can use the standard Monte Carlo simulation 



16 



where Y tjS (s > t) is the stochastic flow of X and is given by 

(**,«);■ = 5 J + J t U (Y t , s )';{d k r\s,X s )ds + <9 fe a l ( S ,X s ) • dW s } (5.23) 

The second term of is already given as an explicit function of x%. 
2nd order 

We can proceed to higher orders in similar fashion. For the second order, the contribution 
to from the first line of can be calculated in the same way as the decoupled case. 
Let us consider non-trivial remaining terms. The contribution from diV^(t,x)/j, l ^°\t,x), 
for example, can be calculated as 



= l{T!>t}E 



1 {r 1 <T}At ( ° ) (n,^r 1 )^-(l{^>r 1 }E 



1 {t2<t}G^(t 2 ,X T2 ) T Tl ) 



l{n<r 2 <T}At (0) (n,X Tl )(y Tl , T2 )^ J G«(r 2 ,X T 



Ft 



where 



(0) 



(s,X s ) = ^eft x " du v m (s,X s ) 



A 



(5.24) 
(5.25) 



Note that the partial derivative of x in diV^ (ti, X Tl ) should be recognized as the shift of 
X at the time of t±, which leads to the first order stochastic flow Y in the above expression. 

Next, let us consider the contribution from dijV^(t, x)(a l ■ r]^ ^)(t,x). As is the pre- 
vious example, it is calculated as 



l{n>t} E 



l {n <T}(^-^°>)(Tl,X n ) 



d 2 



dx l T1 dx 3 T1 



1 {r 2 >r 1 } IE 1 {r 2 <T}G^ > (T2, X T2 ) 



T 



= l{n>t}E [l {T1<T2<T} (a^))(r 1 ,X r j{(r Tl , T J^ fc G'«(r 2 ,X T2 ) 



+ (y T1 , r2 )f(y ri , T2 )^ H G«(r 2 ,X r2 )}|^] (5.26) 



where a 1 ■ is defined similarly as G^ 1 ). Note that the second order stochastic flow 
(Tt,s)ij is defined, for (u > t), as 



(5.27) 



and is given by 



(T 



J^Tt^dir^s, X s )ds + d,<7 fc (s,X s ) • dW a } 

(y M )^(y t , s )™{a /m r fc ( s ,x s )d s + a /m a fc ( s ,x s ) • diy s } (5.28) 



The remaining contributions to as well as can be calculated by the same tech- 
nique. 
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Although tedious calculation is needed, we can proceed to an arbitrary higher order in the 
same fashion. Note that, also in fully-coupled cases, new particles required in simulation 
are all derived as stochastic flows of X and hence the initial values at their creations are 
known beforehand. 

6 Conclusions and Discussions 

In this paper, we have developed an efficient Monte Carlo scheme with an interacting 
particle representation. It allows straightforward numerical implementation to solve fully 
non-linear decoupled as well as coupled FBSDEs at each order of perturbative expansion. 
The appearance of unknown backward components in the expressions of higher order ap- 
proximations is solved by introducing an appropriate particle interpretation. Although a 
couple of new particles are created at random interaction times, their initial values are 
known beforehand. This is due to their properties as the stochastic flows of the under- 
lying sate, which is the crucial point to make straightforward Monte Carlo simulation 
possible. The proposed method can be applied to semi-linear problems, such as American 
and Bermudan options, Credit Value Adjustment (CVA), and even fully non-linear issues, 
such as the optimal portfolio problems in incomplete and/or constrained markets, feed- 
backs from large investors, and also the analysis of various risk measures. It looks also 
interesting to use the current method to study higher order FBSDEs, where the higher 
order Malliavin derivatives exist in the non-linear driver, such as f(t,X t ,Vt,T> t V,'D^V). 
It can be done straightforwardly by introducing higher order stochastic flows. 

Acknowledgment: The authors thank Seisho Sato of the Institute of Statistical Mathe- 
matics (ISM) for the helpful discussions about the branching diffusion method. 
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